Monocyte Subtypes

Setup

# Import functions
root = "./"
source(file.path(root,"general_functions.R"))
import_parameters(params)

## **** __Utilized Cores__ **** = 2$subsetGenes
## [1] "protein_coding"
## 
## $subsetCells
## [1] 500
## 
## $resolution
## [1] 0.6
## 
## $resultsPath
## [1] "./Results"
## 
## $nCores
## [1] 2
## 
## $perplexity
## [1] 30

# load(file.path(resultsPath, "scRNAseq_results.RData"))
resultsPath <- "Results/subsetGenes-protein_coding__subsetCells-F__Resolution-0.2__perplexity-40__nCores-4"
load(file.path(resultsPath,"scRNAseq_results.RData"))

Load Libraries

library(Seurat)
library(cowplot)
library(ggplot2)
library(dplyr) 
library(data.table) 
library(readxl) 
library(reshape2)
library(ggrepel)
library(plotly)

library(GeneOverlap) # BiocManager::install("GeneOverlap")
library(monocle) # BiocManager::install("monocle") 
# BiocManager::install(c('DelayedArray', 'DelayedMatrixStats', 'org.Hs.eg.db', 'org.Mm.eg.db'))
library(org.Hs.eg.db)
library(garnett) # devtools::install_github("cole-trapnell-lab/garnett")
library(enrichR) #BiocManager::install("enrichR")

Identify Cell Types with Garnett + Monocle

Pre-trained PBMC Classifer

# Convert Seurat objt to CDS object 
mDAT <- monocle::importCDS(DAT,  import_all = T) 
# generate size factors for normalization later
mDAT <-  DESeq2::estimateSizeFactors(mDAT)
# Get pre-trained PBMC classifer 
load(file.path(root, "Data/Garnett/hsPBMC")) # Download from: https://cole-trapnell-lab.github.io/garnett/classifiers/hsPBMC  

mDAT <- garnett::classify_cells(mDAT, hsPBMC,
                           db = org.Hs.eg.db,
                           cluster_extend = TRUE,
                           cds_gene_id_type = "SYMBOL")

## Warning in doTryCatch(return(expr), name, parentenv, handler): The
## following genes used in the classifier are not present in the input CDS.
## Interpret with caution. ENSG00000174059The following genes used in the
## classifier are not present in the input CDS. Interpret with caution.
## ENSG00000007062The following genes used in the classifier are not present
## in the input CDS. Interpret with caution. ENSG00000157404The following
## genes used in the classifier are not present in the input CDS. Interpret
## with caution. ENSG00000185291

table(pData(mDAT)$cell_type)

## 
##         B cells           CD34+     CD4 T cells     CD8 T cells 
##              13              92               3               1 
## Dendritic cells       Monocytes        NK cells         T cells 
##             280           15025             801              16 
##         Unknown 
##            2913

table(pData(mDAT)$cluster_ext_type) 

## 
##         B cells           CD34+     CD4 T cells Dendritic cells 
##              12              17               3             165 
##       Monocytes        NK cells         T cells         Unknown 
##           17165             701              15            1066

# Get feature genes for each cell type
feature_genes <- garnett::get_feature_genes(classifier = hsPBMC,
                                   node = "root",
                                   db = org.Hs.eg.db,
                                   convert_ids = F)
head(feature_genes) 

##                      B cells       CD34+ Dendritic cells     Monocytes
## (Intercept)     -0.884989922  0.47694427    -0.124629734 -0.0301907756
## ENSG00000196154 -0.005499815 -0.01088859    -0.007978714  0.0056365451
## ENSG00000019582  0.030146634 -0.01959215     0.015232191  0.0001883507
## ENSG00000156738  2.193353690 -0.66191294    -1.104899464 -0.7393680112
## ENSG00000177455  2.056641938 -0.71851810    -1.941859710 -1.0448682467
## ENSG00000105369  1.800296145 -0.53835614     0.063380804 -0.1973728377
##                     NK cells      T cells     Unknown
## (Intercept)     -0.911657164  0.524757370 0.949765956
## ENSG00000196154 -0.001627233 -0.005835858 0.026193667
## ENSG00000019582 -0.005880096 -0.023473024 0.003378099
## ENSG00000156738 -0.096566144 -0.602417257 1.011810124
## ENSG00000177455  0.297279508  0.177790907 1.173533704
## ENSG00000105369 -0.618567693 -0.963141845 0.453761563

# Convert back to Seurat object & save results
DAT <- Seurat::AddMetaData(DAT, pData(mDAT)[c("garnett_cluster","cell_type","cluster_ext_type")])
write.csv(DAT@meta.data, file.path(resultsPath, "garnett_results.csv"), quote = F,row.names = F)

Biomarker Expression

markerList <- c("CD14", "FCGR3A")
markerData <- DAT@scale.data[row.names(DAT@scale.data) %in% markerList,] %>% t() %>% data.frame()  
## Efficiently merge using data.table
dt1 <- data.table(markerData, keep.rownames = "Cell", key = "Cell") %>% copy()
dt2 <- data.table(DAT@meta.data[c("cell_type","post_clustering")], keep.rownames = "Cell", key = "Cell") %>% copy()
markerDT <- dt1[dt2]
 
ggplot(data = markerDT, aes(x=CD14, y=FCGR3A) ) + 
  stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="purple", bins = 100, size=.1) +
  geom_point(aes(color=as.factor(cell_type)), shape=16, stroke=0, size=2, alpha=.1) + 
  guides(colour = guide_legend(title="cell_type",override.aes = list(alpha = 1))) 

tSNE Metadata Plots

tSNE_metadata_plot <- function(var, labSize=12){ 
  # Metadata plot  
  p1 <- Seurat::TSNEPlot(DAT, do.return = T,  do.label = T,  group.by = var,label.size = labSize,
                 plot.title=paste(var), vector.friendly=T) +
    theme(legend.position = "top") + 
    scale_color_brewer( palette = "Dark2",  aesthetics = aes(alpha=.5)) 
     
  # t-SNE clusters plot
  p2 <- Seurat::TSNEPlot(DAT, do.return = T, do.label = T,label.size = labSize,
                 plot.title=paste("Unsupervised Clusters"), vector.friendly=T)  +
    theme(legend.position = "top") + 
    scale_color_brewer( palette = "Set1", aesthetics = aes(alpha=.5))  
  print(cowplot::plot_grid(p1,p2))
}   

Garnett cell_type

tSNE_metadata_plot("cell_type")

Garnett cluster_ext_type

tSNE_metadata_plot("cluster_ext_type")

mut

tSNE_metadata_plot("mut")

tSNE FeaturePlots

Monocyte Markers

FeaturePlot(DAT,features.plot =  c("CD14","FCGR3A"),  
            cols.use = c("grey","red","blue","purple"),
            dark.theme = T, nCol = 2, overlay = T, no.legend = F) 

PD Genes

FeaturePlot(DAT,features.plot = c("LRRK2", "GBA"),  
            cols.use = c("grey","purple","cyan","blue"),
            dark.theme = T, nCol = 2, overlay = T, no.legend = F) 

MS4A4A & MS4A6A

FeaturePlot(DAT, features.plot =  c("MS4A4A","MS4A6A"), 
            cols.use = c("grey","red","green","blue"),
            dark.theme = T, nCol = 2, overlay = T, no.legend = F)

AD/PD-related Gene Expression

Q: Are AD/PD-related gene expression more prevalent in classical or intermediate monocyte subtypes?

# Combine AD gene lists
curatedGenes <- read_excel("Data/curated_AD-PDgene_lists.xlsx")
AD_related_genes <- read_excel(file.path(root, "Data/AD-related_genes.xlsx"))
ADgenes <- rbind(AD_related_genes, data.frame(Gene=curatedGenes$AD_panel, Category="AD") ) %>% unique()
# Combine PD gene lists
PDgenes <- data.frame(Gene=c("LRRK2","GBA",
                             curatedGenes$gwas_Nearest_gene_Nalls_2019, 
                             curatedGenes$QTL_Nominated_genes_Nalls_2019,
                             curatedGenes$TWAS_Garrett
                             ), Category="PD") %>% unique()


MonocyteSubtype.markers <- read.csv(file.path(root, resultsPath, "MonocyteSubtype.markers.csv"), row.names = 1)
MonocyteSubtype.markers_sig <- subset(MonocyteSubtype.markers, p_val_adj <= 0.05) 
createDT(MonocyteSubtype.markers, "Full list of DEGs: Cluster 0 vs. Cluster 1")

Test Gene Overlap

Functions

Gene Overlap Functions

• “The GeneOverlap class formulates the problem as testing whether two variables are independent, which can be represented as a contingency table, and then uses Fisher’s exact test to find the statistical significance.”
• Low p-value means the overlap is highly significant.
• Documentation

report_overlap <- function(DAT, list1, list2){ 
  genomeSize <- dim(DAT@raw.data)[1]
  go.obj <- newGeneOverlap(listA = list1, listB = list2, genome.size = genomeSize )
  go.obj <- testGeneOverlap(go.obj)
  print(go.obj) 
  
  
  overlappingGenes <- getIntersection(go.obj)
  percent_of_targetGenes <- length(overlappingGenes) / length(list2)*100
  percent_of_DEGs <- length(overlappingGenes) / length(list1)*100
  targetGenes_DEGs <- list2[overlappingGenes %in% list2]
  
  cat("\n",round(percent_of_targetGenes, 2),"% of the provided genes (",length(overlappingGenes),"/",length(list2),") are differentially expressed between Canonical vs. Intermediate monocyte subtypes.")
   
  cat("\n",round(percent_of_DEGs, 2),"% of the DEGs between Canonical vs. Intermediate  monocyte subtypes (",
      length(overlappingGenes),"/",length(list1),") are in the provided gene list.")
  
  cat("\n-------------------------------------------------------",
      "\n\n********** Enrichment p-value =",go.obj@pval,"**********\n\n")
  return(overlappingGenes) 
} 

Overlap Plot Functions

# PLOT FUNCTIONS
overlap_heatmap <- function(clustDAT, geneList, title="Overlapping Genes:\nGene List vs. DGE Genes"){ 
  markerDF  <- get_markerDT(clustDAT, markerList = geneList, rawData = T)
  markerMatrix <- reshape2::acast(markerDF, Gene~post_clustering, value.var="Expression",
                                  fun.aggregate = mean, drop = F, fill = 0)
  # Heatmap.2 
  my_palette <- colorRampPalette(c("purple", "black", "cyan"))(n = 1000)
  hmap <- gplots::heatmap.2(markerMatrix, xlab = "Cluster", dendrogram = "row",
                    col = my_palette, tracecol = "gray", srtCol = 0, offsetCol=1.5, vline=T, 
                    trace='none', key.title=NA, key.ylab = "Expression", colsep=T, sepwidth = 0.01, 
                    main = title) 
}


make_fractionDF <- function (MonocyteSubtype.markers_sig, geneList, allGenes=F){
  MonocyteSubtype.markers_sig$Gene <- row.names(MonocyteSubtype.markers_sig)
 if(allGenes==F){
   geneDF <- subset(MonocyteSubtype.markers_sig, Gene%in%geneList)    
 } else{geneDF <- MonocyteSubtype.markers_sig}
  pct_df <- melt(geneDF, id.vars = c("Gene","avg_logFC"), measure.vars = c("pct.2", "pct.1"), 
                 variable.name = "Cluster", value.name = "FractionCells") 
  pct_df$Cluster <- ifelse(pct_df$Cluster=="pct.1", "1", "0") 
  return(pct_df)
}

cellFractions_plot <- function (MonocyteSubtype.markers_sig, geneList, title="", allGenes=F){ 
  pct_df <- make_fractionDF(MonocyteSubtype.markers_sig, geneList, allGenes) 
  # Fraction Cells
  cfp <- ggplot(data=pct_df, aes(x=Gene, y=FractionCells, fill=Cluster)) + geom_col(position="dodge") + 
    labs(title = title, y="Fraction of Cells", x="Gene") +  
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
  if(allGenes==T){
    cfp <- cfp +  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8))
    print(ggplotly(cfp))
  }else{print(cfp)}
}

logFC_plot <- function(MonocyteSubtype.markers_sig, geneList, title="", allGenes=F){
  pct_df <- make_fractionDF(MonocyteSubtype.markers_sig, geneList, allGenes) 
  # LogFC
  lfcp <- ggplot(data=pct_df, aes(x=Gene, y=avg_logFC, fill=Gene)) + geom_col(position="dodge") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
     labs(title = title, y="avg_logFC", x="Gene")
   if(allGenes==T){
    lfcp <- lfcp +  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8))
    print(ggplotly(lfcp))
  }else{print(lfcp)} 
} 

enrichment_test_and_plots <- function(DAT, DEG_df, geneList, disease=""){
  overlappingGenes <- report_overlap(clustDAT, 
                                      list1=row.names(DEG_df), 
                                      list2=geneList)
  ## All DEGs
  # cellFractions_plot(DEG_df, geneList,
  #                      title= paste("Fraction of Cells Expressing\nAll",disease,"Genes"), allGenes = T)
  # logFC_plot(DEG_df, geneList,  title= "LogFC\nAll PD Genes", allGenes = T)
  # Overlapping Genes
  if(length(overlappingGenes)>0){
    FeaturePlot(DAT, features.plot = overlappingGenes, dark.theme = T) 
    cellFractions_plot(DEG_df, overlappingGenes,
                       title= paste("Fraction of Cells Expressing\n",disease,"Genes Overlapping with DEGs") )
    logFC_plot(DEG_df, overlappingGenes, title= paste("LogFC\n",disease,"Genes Overlapping with DEGs") )
    
  
    if(length(overlappingGenes)>1){
      overlap_heatmap(clustDAT, geneList=overlappingGenes, 
                      title=paste("Overlapping Genes:\n",disease,"Genes vs. DGE Genes") )
    } 
  } else {cat("There were no overlapping genes between the DGE list and",disease,"genes.")} 
  return(overlappingGenes)
}


clustDAT <- SubsetData(DAT, subset.name = "post_clustering", accept.value = c(0,1), do.scale = F) 
purpleScale <-  c("grey","purple","blueviolet","magenta")
heatScale <-  c("yellow","grey","red") 

PD Genes vs. DGE Overlap

overlappingGenes_PD <- enrichment_test_and_plots(clustDAT, MonocyteSubtype.markers_sig, PDgenes$Gene)  

## Detailed information about this GeneOverlap object:
## listA size=205, e.g. FCGR3A RHOC CDKN1C
## listB size=134, e.g. LRRK2 GBA HRNR
## Intersection size=1, e.g. GRN
## Union size=338, e.g. FCGR3A RHOC CDKN1C
## Genome size=24914
## # Contingency Table:
##       notA inA
## notB 24576 204
## inB    133   1
## Overlapping p-value=0.67
## Odds ratio=0.9
## Overlap tested using Fisher's exact test (alternative=greater)
## Jaccard Index=0.0
## 
##  0.74 % of the provided genes ( 1 / 135 ) are differentially expressed between Canonical vs. Intermediate monocyte subtypes.
##  0.49 % of the DEGs between Canonical vs. Intermediate  monocyte subtypes ( 1 / 205 ) are in the provided gene list.
## ------------------------------------------------------- 
## 
## ********** Enrichment p-value = 0.6704846 **********

AD Genes vs. DGE Overlap

overlappingGenes_PD <- enrichment_test_and_plots(clustDAT, MonocyteSubtype.markers_sig, ADgenes$Gene)  

## Detailed information about this GeneOverlap object:
## listA size=205, e.g. FCGR3A RHOC CDKN1C
## listB size=76, e.g. SPI1 CD33 PILRA
## Intersection size=10, e.g. IFITM3 IFITM2 IFITM1
## Union size=271, e.g. FCGR3A RHOC CDKN1C
## Genome size=24914
## # Contingency Table:
##       notA inA
## notB 24643 195
## inB     66  10
## Overlapping p-value=6.8e-10
## Odds ratio=19.1
## Overlap tested using Fisher's exact test (alternative=greater)
## Jaccard Index=0.0
## 
##  12.82 % of the provided genes ( 10 / 78 ) are differentially expressed between Canonical vs. Intermediate monocyte subtypes.
##  4.88 % of the DEGs between Canonical vs. Intermediate  monocyte subtypes ( 10 / 205 ) are in the provided gene list.
## ------------------------------------------------------- 
## 
## ********** Enrichment p-value = 6.795736e-10 **********

Identify Cluster-specific Biomarkers

clust0.markers <- FindMarkers(DAT, ident.1=0, min.pct = 0.25, only.pos = F, test.use = "wilcox")
clust1.markers <- FindMarkers(DAT, ident.1=1, min.pct = 0.25, only.pos = F, test.use = "wilcox")

clust0.uniqueMarkers <- clust0.markers[!(row.names(clust0.markers) %in% row.names(clust1.markers)),] %>%
  subset(p_val_adj<=0.05)
clust1.uniqueMarkers <- clust1.markers[!(row.names(clust1.markers) %in% row.names(clust0.markers)),] %>%
  subset(p_val_adj<=0.05)

difference <- abs( length(row.names(clust0.uniqueMarkers)) - length(row.names(clust1.uniqueMarkers) ) )
uniqueMarkers <- data.frame(Cluster0_markers=c(row.names(clust0.uniqueMarkers), rep("",difference) ),
                            Cluster1_markers=row.names(clust1.uniqueMarkers))
write.csv(uniqueMarkers,
          file.path(resultsPath,"unique_cluster_markers.csv"), 
          quote = F, row.names = F)

DGE Across Clusters

DGE Across Clusters: dx

# clustDAT@meta.data$dx %>% unique()
DEGs <- runDGE(clustDAT, meta_var = "dx", group1 = "PD", group2 = "control")

DGE Across Clusters: mut

# clustDAT@meta.data$dx %>% unique()
DEGs <- runDGE(clustDAT, meta_var = "mut", group1 = "PD", group2 = "GBA")

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_text_repel).

DGE Within Clusters

DGE Within Clusters: dx

DGE_within_clusters(clustDAT,  meta_var = "dx", group1 = "PD", group2 = "control", clusterList = c(0,1))

Cluster 0: PD vs. control

Cluster 1: PD vs. control

DGE Within Clusters: mut

DGE_within_clusters(DAT, meta_var = "mut", group1 = "PD", group2 = "GBA", clusterList = c(0,1))

Cluster 0: PD vs. GBA

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_text_repel).

Cluster 1: PD vs. GBA

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_text_repel).

DEG Enrichment w/ enrichR

enrichr_dbs <- c("KEGG_2018", "Reactome_2016",
                 "GO_Biological_Process_2018", "GO_Molecular_Function_2018", "GO_Cellular_Component_2018", 
                 "Rare_Diseases_AutoRIF_ARCHS4_Predictions", "Human_Gene_Atlas")
# createDT(enrichR::listEnrichrDbs(), "Enrichr Databases")
 
geneList <- data.frame(Gene=row.names(MonocyteSubtype.markers), 
     Weight=scales::rescale(length(MonocyteSubtype.markers$p_val_adj):1))

results <- enrichr(genes = geneList, databases = enrichr_dbs ) 

Uploading data to Enrichr… Done. Querying KEGG_2018… Done. Querying Reactome_2016… Done. Querying GO_Biological_Process_2018… Done. Querying GO_Molecular_Function_2018… Done. Querying GO_Cellular_Component_2018… Done. Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions… Done. Querying Human_Gene_Atlas… Done. Parsing results… Done.

for (db in enrichr_dbs){
  cat('\n')
  cat("##",db,"\n")  
  # res <- subset(results[[db]], Adjusted.P.value<=0.05)
  createDT_html(results[[db]], paste("Enrichr Results:",db))
  cat('\n')
}  

KEGG_2018

Reactome_2016

GO_Biological_Process_2018

GO_Molecular_Function_2018

GO_Cellular_Component_2018

Rare_Diseases_AutoRIF_ARCHS4_Predictions

Human_Gene_Atlas